Symbiotic Self-Organized Criticality and Intermittent Turbulence 
in an MHD Current Sheet with a Threshold Instability 
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We report numerical evidence of a self-organized criticality (SOC) and intermittent turbulence 
(IT) symbiosis in a resistive magnetohydrodynamic (MHD) current sheet model that includes a 
local hysteretic switch to capture plasma physical processes outside of MHD that are described 
in the model as current-dependent resistivity. Results from numerical simulations show scale-free 
avalanches of magnetic energy dissipation characteristic of SOC, as well as multi-scaling in the 
velocity field numerically indistinguishable from certain hierarchical turbulence theories. We suggest 
that SOC and IT may be complementary descriptions of dynamical states realized by driven current 
sheets - which occur ubiquitously in astrophysical and space plasmas. 
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Most theoretical studies of self-organized criticality 
(SOC) focus on cellular models such as the paradigmatic 
BTW sandpile [2, Q. On the other hand, many ex- 
amples of SOC-like phenomena occur in systems whose 
canonical descriptions invoke continuum equations such 
as the Navier-Stokes or magnetohydrodynamic (MHD) 
equations, whose solutions can exhibit some properties 
of intermittent turbulence (IT) . It has been argued that 
complementarity between SOC and IT can be realized 
in avalanching systems 0, 0, H, H, 0, @, • The pro- 
posal that SOC and IT could be distinguished by ana- 
lyzing waiting times between bursts [l(| has turned out to 
be false: Once a finite observational threshold (unavoid- 
able in any physical measurement) is introduced, even 
the ordinary BTW sandpile exhibits scale free waiting 
time statistics [7j. Simultaneous signatures of SOC and 
IT have been observed in physical systems such as, e.g., 
the solar corona through an analysis of an extended set 
of high resolution images provided by the SOHO space- 
craft [8(. Earth's magnetosphe re p roduces power-law dis- 
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as well as related 
Avalanches 



intermittent plasma turbulence 
with signatures of critical behavior have also been de- 
tected in numerical fluid models 17, 18[ but the poten- 



tial complementarity between SOC and IT has not been 
rigorously proved. 

Here, we analyze a current sheet model 



18, UM that 



supports a magnetic field reversal in an MHD plasma. 
The model contains measurable plasma parameters cou- 
pled through the full system of resistive MHD equations, 
with boundary conditions that allow a persistent cur- 
rent sheet to form. To describe local breakdown, where 
standard MHD equations fail and the true plasma equa- 
tions must be invoked, we adapt the idea of E. Lu [20|. 
A current dependent threshold instability, taking the 
form of an hysteretic switch controlling plasma resistiv- 
fl9j |. is used. With a uniform drive, the model 
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ing phases. We observe simultaneous signatures of SOC 
and IT in the dynamics of the model during the unload- 
ing phase. We limit our discussion to two of these: (1) 
We show that the model exhibits scale-free avalanching, 
consistent with SOC, of magnetic energy into a central 
reconnection zone where the field is dissipated, and (2) 
we provide evidence that this avalanching drives fluc- 
tuations in the velocity field that exhibit multi-scaling 
indistinguishable from certain hierarchical models of IT 
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repeatedly cycles through large scale loading and unload- 



Klimas ct al. [2004, 2005] have shown that the presence 
of the hysteretic switch in the model leads to the propaga- 
tion of micro-scale current waves that induce concurrent 
waves of resistivity, thus leading to local unfreezing of 
the magnetic field to form field energy avalanches with 
scale-free statistics. Here we provide evidence that lo- 
calized, intermittent J x B acceleration associated with 
these current waves is directly responsible for the multi- 
scaling velocity field. In this current-sheet model, the 
SOC and IT dynamics are intimately associated in this 
way. We have found it impossible to have one without 
the other. 

The 2d numerical model discussed here was originally 
motivated by observations of Earth's magnetotail, which 
is a driver for high-latitude geomagnetic activity [23j |. 
The key plasma structure that controls geomagnetic sub- 
storms and is represented by our model is the large-scale 
cross-tail current sheet. The model consists of the full 
compressible, resistive MHD system, including a poly- 
tropic energy equation (see eq.(l) - (4) in [181 ] for a com- 
plete definition). In addition the local resistivity, D, 
obeys the following equations: 

Q = { D min , \J\ < PJ C ; D max ,\J\>J c } (1) 
d t D= (Q(\J\)-D)/r , (2) 

which were adapted from Lu [2fJ. The function Q takes 
the value D max ks 1 wherever the local current density J 
exceeds a critical value J c and returns to D m i n <C D max 
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when the current density falls below (3J C , with j3 < 1. 
Hence, the switch in Q is hysteretic. The transition from 
Dmin t° D max represents the excitation and saturation 
of a kinetic current-driven instability over a time interval 
that is below the resolution of MHD and, hence, enters as 
an instantaneous transition. The transition from D max 
to D m i n represents the subsequent quenching of the in- 
stability. The resistivity, D, is assumed to grow or decay 
with a single time-scale r that is slow compared to the 
simulation time step. The values of the model parame- 
ters used here can be found in 
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We have obtained numerical solutions of the current 
sheet model on a 400 x 400 grid with D m i n <C 1 such 
that wherever D = D m i n the evolution is indistinguish- 
able from ideal MHD over the observed time scales. A 
configuration of the model in a snapshot during an un- 
loading phase is illustrated in Fig. 1. Plasma inflowing 
at the upper and lower boundaries carries magnetic flux 
with it, thus increasing the strength of the magnetic field 
reversal (and hence the electric currents in the current 
sheet). Eventually the current reaches J = J c somewhere 
- not necessarily at z = - which starts an avalanche of 
magnetic flux transport toward the central region where 
flux is annihilated and converted to kinetic and thermal 
energy. This conversion process also drives plasma out 
of the region through the open boundary at the right. A 
small portion (10 -3 - 10~ 2 ) of the input magnetic energy 
is carried out through this boundary as well. 
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FIG. 1: Numerical simulation of a current sheet and its quali- 
tative behavior (see pl| for details). The model is driven by a 
steady, uniform plasma inflow at the top and bottom bound- 
aries as shown by arrows. Fhe left boundary is closed, the 
right is open. Plasma energized through annihilation in the 
magnetic field reversal leaves the region at the right. Inset: 
(Top) A snapshot of regions where the diffusive Poynting flux 
exceeds a threshold, used to define avalanches. (Bottom) At 
the same time, the corresponding velocity field (v x ) in the sys- 
tem. Note that while the magnetic field lines appear smooth 
at this scale, the plasma velocity field is highly intermittent. 

Eventually, the simulated plasma reaches a statistically 
stationary state in which the rates of magnetic energy 



and plasma mass flowing into the region are balanced, 
over long time scales, by the field annihilation rate and 
the outflow at the open boundary. After about 10 3 Alfven 
traversal times, this state takes the form of large scale 
global cycles consisting of long laminar periods during 
which plasma is loaded into the system but no active 
grid sites (with Q = D max ) are generated, followed by 
highly erratic unloading periods during which the mag- 
netic field undergoes local transitions between frozen and 
unfrozen states analogous to stick-slip behavior of SOC 
models 24, 25, [26|, [27|. Fig. 2 illustrates the initiation 
of one of these unloading periods. Neglecting a small 
convective contribution, the magnetic energy transported 
when a site becomes active is given by the local diffusive 
Poynting flux Sa(x,z) = [cjAa:)r\J x B [l8j, in which 
rj = 2irD/c 2 is the anomalous resistivity, and c is the 
speed of light. Fig. 2 shows the initial expansion of a 
Poynting flux avalanche in back of an outward propagat- 
ing current wave as well as the associated transition from 
laminar to turbulent plasma flow at and in the interior 
of the expanding wave. 

To observe avalanches, we used an automated tech- 
nique for detecting and tracing regions having grid sites 
with Sd above a certain threshold. In analogy with 
avalanches in 2d sandpiles, we treat the events as 2+1 di- 
mensional spatiotemporal objects. Avalanches were iden- 
tified by applying a floating activity threshold Sth(t) = 
(Sd) + k ■ a adjusted to the average value (Sd) and the 
standard deviation a of the Poynting flux at every time 
step [3^]. The time evolution of avalanches was ob- 
tained by checking the intersections of spatial regions 
above Sth(t) in consecutive pairs of Sd(x,z) snapshots. 
Each avalanche was characterized by its lifetime T and 
its total Poynting flux, E, obtained from the integra- 
tion of Sd over its spatiotemporal domain - grid sites 
with Sd(x,z,t) > Sth(t) taking part in the avalanche. 
The linear dimensions l x and l z of avalanches were esti- 
mated by determining standard deviations of the x and 
z coordinates over all the grid sites involved in each 
avalanche (equally weighted). In addition the geometric 
mean l xz = (Ixlz) 1 ^ 2 , as well as the total area s, rep- 
resenting the total number of distinct pixels involved in 
the avalanche, were estimated. The statistics reported 
here were obtained using k = 3.0 and have also been 
reproduced in the range k = 1.5 — 4.0. 

The probability distribution for lengths, time, area, 
and energy of avalanches all obey scale free statistics. 
The first group of critical exponents was estimated based 
on analyses of probability distributions p(T, s max ) and 
p(E,s max ) constructed from subsets with s < s max , in 
which s max is defined to be the maximum area (number 
of pixels) of events included in the subset used to make 
the histogram. The normalized probability distributions 
were studied using the scaling ansatz 
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fx(X/X c ), X c ~ 8 >_ 
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FIG. 2: The neighborhood of the initial site of instability 
shortly after the initiation of an unloading event, (a) Poynt- 
ing flux due to slipping magnetic flux in the interior of the 
outward expanding current wave, (b) Transition from lami- 
nar to turbulent velocity field due to J x B acceleration at 
the current wave as it propagates through the magnetic field. 



where X E {E, T} and fx are scaling functions that are 
approximately constant for X < X c and drop rapidly for 
X > X c . Assuming Eq. [31 we have plotted the distri- 
butions in the rescaled coordinates (X / s^ ax , p(X)X Tx ) 
and identified the combination of te, tt, Xe and At ex- 
ponents that provides the best data collapse (Fig. 3). The 
resulting values t e = 1.48±0.02 and t t = 1.95±0.03 co- 
incide with those reported earlier in [1 81 ] - The exponents 
X E = 1.47±0.03 and A T = 0.68±0.04 are consistent with 
the regression analyses for the expected values of energy 
and lifetime for avalanches with a given size s {E) s and 
(T) s , within statistical error. These values also preserve 
the scaling relation Xt{tt — 1) = Xe(te — 1)- 

The anisotropy of the model leads to different growth 
rates of avalanches in the x and z directions. Hence 
~ (U dx and (s) iz ~ (l z ) d * with 4 = 1.40±0.03and 
d z = 3.11 ± 0.06. The geometric mean l xz is related to s 
through another scaling relation, (s); xz ~ (l X z) d:c * , with 
d X z — 1-97 ± 0.05 indicating that avalanches are com- 
pact. The exponent values obtained are consistent with 
d-x 1 +d~ 1 = 2d~}. We also studied E and T as functions 
of l x and l z and found that {E) lx ~ 1%*, (T) /<e - l»<° , 
(E)i z ~ 1% Z and (T)i z ~ l v z z with scaling exponents 
fi x = 1.87 ± 0.04, v x = 1.13 ± 0.01, [i z = 3.50 ± 0.04 
and v z = 1.97 ±0.06 (see the insets in Fig. 3). The ratios 
M^/Mz and v x /v z are close to d x /d z as expected. 

All these results indicate that with respect to bursts 
of energy dissipation above background, the system op- 
erates at or near a SOC state. Based on the values of 
the distribution exponents te and tt one can conjecture 
that the model operates near the mean-field limit. This is 
not typical for non-directed SOC sandpiles whose upper 
critical dimension d u is usually higher than 2. However, 
it is possible that the avalanching dynamics in our model 
can be mapped onto the universality class of non-Abelian 
directed sandpiles with irreversible topplings [28] which 
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FIG. 3: Data collapse using Eq. 3 for avalanches with different 
maximum area s m ax with tt = 1.95, At = 0.68, te = 1.48, 
and Xe = 1.47. Insets: anisotropic scaling of E and T with 
the maximum avalanche extent I in x and z directions. 



exhibit the mean-field exponents starting from d u =2. 

To analyze the velocity field, we have computed a set of 
equal time structure functions defined as S q (l) — {\Svi\ q ), 
where 5vi = (v(r + 1) — v(r)) ■ I / 1 is the increment of the 
velocity v in the direction I (parallel to x or z axes), q is 
the order of the structure function, I = \l\ is the spatial 
displacement, and averaging indicated by (• • ■) is per- 
formed over all positions r and times during an unload- 
ing phase. For many turbulent phenomena, S q (l) ~ Z^ 9 ' 
with defined by the turbulent regime under study. 
To extend the scaling range and improve the accuracy 
of this analysis, we have applied the method of extended 
self similarity (ESS) [2{J by plotting S q (l) versus S 3 (l). 

The resulting structure functions (Fig. 4) exhibit ESS 
over the entire range of scales available. The error bars 
shown are for C(<?)/C(<3) in the x direction; the errors in 
the z direction are about three times smaller. The values 
obtained in both directions are the same up to these er- 
rors. The dependence of C(<?)/C(3) on the order q shows a 
systematic departure from the Kolmogorov law £ = q/2>. 
In principle, it can be fitted by the hierarchical model 
C(q) = (1 - l)q/g + C(l - [1 - 7/C] 9 / 9 ], in which C is 
the codimension of the most singular dissipative struc- 
tures, g and 7 are defined by 6 z ~ i 1 ^ 9 and t e ~ P , 
with t e being the energy transfer time at the smallest in- 
ertial scales £ [3(J. By choosing either g = 4, 7 = 1/2, 
C = 1 (Iroshnikov-Kraichnan theory (IK) j22j) or g = 3, 
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7 = 2/3, C = 2 (She and Leveque (SL) theory [21() one 
can obtain rather accurate fits to the data. However, the 
physical conditions for the turbulence in our model are 
strikingly different from those in any of these hierarchical 
models. The current waves at the leading fronts of energy 
avalanches accelerate the fluid through the J x B force. 
These current structures play the role of energy sources 
rather than energy sinks (as would be the case in classi- 
cal turbulent models). As we have shown, the avalanches 
are scale-free, and thus this driving mechanism appears 
at all scales as opposed to the standard picture of the 
direct turbulent cascade. At this point there is no closed 
theory that could describe the relationship between the 
magnetic energy avalanches and the intermittency in the 
velocity field (e.g. in the form of an "exact law" [3lj]). 
although it is evident that the resulting behavior is dif- 
ferent from the MHD turbulence 13911. 
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FIG. 4: Left: ESS plots of velocity structure functions with 
I parallel to z axis. The inset shows the same functions with 
subtracted average slopes. Right: Dependence of C(<z)/C(3) 
on the order q for our current sheet model (CSM, with error 
bars), hierarchical models of IT mentioned in the text, Muller 
and Biskamp (MB) [H model (g = 3, 7 = 2/3, C = 1), 
Kolmogorov model (K41) 33|, as well as the exponents from 
2d ideal MHD simulations [34l |. 

The interplay between SOC- and turbulence-based 
routes to multiscale complexity in natural systems re- 
mains a subject of extreme interest in the modern phys- 
ical literature. There is a growing body of evidence that 
statistical signatures of these routes coexist in real-world 
processes, including those in the solar corona and Earth's 
magnetosphere. At the same time, there is no clear un- 
derstanding of how such coexistence arises from first prin- 
ciples of fluid dynamics. Our study partially fills this 
gap by providing an example of a realistic (and very 
common, see e.g. | 



361 [37 [) plasma configuration in 



which avalanches of electromagnetic energy stir an oth- 
erwise laminar velocity field to generate a self-consistent 
pattern of multiscale spatiotemporal fluctuations obey- 
ing both SOC and turbulent scaling laws. There is no 
doubt that the resulting behavior is a more general form 



of dynamical complexity than either of the two scenar- 
ios separately. We are looking forward to seeing future 
theoretical studies of this phenomenon. 

The work of A. J. K. was supported by the NASA 
Geospace Sciences program. 



[1 

[2. 
[3. 
[4 
[5 

ir 

[8 
[9 
[10 

[11 

[12 

[13 

[14 
[15. 
[16 
[17; 

[18 

[19 

[20 
[21 

[22 
[23 

[24 

[25 
[26 
[27 
[28 

[29 
[30 
[31 

[32 

[33 
[34 
[35 



[36 

[37 
[38 



Bak et al., Phys. Rev. Lett. 59, 381 (1987). 
Bak et al., Phys. Rev. A 38, 364 (1988). 
Bak and M. Paczuski, Physica A 348, 277 (2005). 
Chen et al., Physica A 340, 566 (2003). 
Chen and C. Jayaprakash, Physica A 340, 566 (2004). 
R. Sreenivasan et al., Physica A 340, 574 (2004). 
Paczuski et al., Phys. Rev. Lett. 95, 181102 (2005). 
M. Uritsky et al., e-print astro-ph/0610130 (2006). 
Chang, Phys. Plasmas 6, 4137 (1999). 
Boffetta et al., Phys. Rev. Lett. 83, 4662 (1999). 
M. Uritsky et al., J. Geophys. Res. -Space Physics 107, 
1426 (2002). 

V. M. Uritsky et al., Geophys. Res. Lett. 30, 1813 (2003). 
V. M. Uritsky et al., Geophys. Res. Lett. 33, L08102 

(2006) . 

V. Angelopoulos et al., Phys. Plasmas 6, 4161 (1999). 

A. T. Y. Lui et al, Annales Geophysicae 24, 2005 (2006). 

Z. Voros et al., Space Sci. Rev. 122, 301 (2006). 

P. Dmitruk and D. O. Gomez, Astrophys. J. 484, L83 

(1997). 

A. J. Klimas et al., J. Geophys. Res.-Space Physics 109, 
1426 (2004). 

A. J. Klimas et al., Geophys. Res. Lett. 32, L14108 
(2005). 

E. T. Lu, Phys. Rev. Lett. 74, 2511 (1995). 

Z.-S. She and E. Leveque, Phys. Rev. Lett. 72, 336 

(1994). 

R. Grauer et al., Phys. Lett. A 195, 335 (1994). 

D. N. Baker et al., J. Geophys. Res.-Space Physics 104, 

14601 (1999). 

P. Bak and C. Tang, J. Geoph. Res.-Solid Earth and 
Planets 94, 15635 (1989). 

Z. Olami et al., Phys. Rev. Lett. 68, 1244 (1992). 
M. Paczuski et al., Phys. Rev. E 53, 414 (1996). 
M. Gedalin et al., Phys. Rev. E 72, 037103 (2005). 

D. Hughes and M. Paczuski, Phys. Rev. Lett. 88, 054302 
(2002). 

R. Benzi et al., Phys. Rev. E 48, R29 (1993). 

H. Politano and A. Pouquet, Phys. Rev. E 52, 636 (1995). 

C. Connaughton et al., Phys. Rev. Lett. 98, 080601 

(2007) . 

W.-C. Muller and D. Biskamp, Phys. Rev. Lett. 84, 475 
(2000). 

A. Kolmogorov, Dokl. Akad. Nauk SSSR 31, 538 (1941). 
H. Politano et al., Europhys. Lett. 43, 516 (1998). 

E. R. Priest and T. G. Forbes, Magnetic Reconnection: 
MHD Theory and Applications (Cambridge Univ. Press, 
2000). 

M. J. Aschwanden, Physics of the Solar Corona 
(Springer, 2005). 

S. C. Cowley et al., Phys. Reports 283, 227 (1997). 
The floating threshold allows improved statistics consid- 
ering the time variation during unloading cycles 



5 

[39] This is also evident from the 3d version of our model 
which will be described in a separate paper. 



